Sie dürfen für die Lösung der Übung zusammenarbeiten, so lange sich ihre Zusammenarbeit auf algorithmische Fragestellungen beschränkt. Bei kopiertem Code oder Text (von Mitstudierenden oder dem Internet) werden alle Lösungen der beteiligten Parteien mit 0 Punkten bewertet. Dazu werden alle Lösungen manuell und automatisiert auf Kopien untersucht.
Einzureichen (online auf mlhub submitten) bis Montag, 2. April 2018, 24:00 Uhr
Ziel dieser Übung ist, den Preis einer Immobilie mit Hilfe von linearen Regressionen aus verschiedenen Attributen abzuschätzen.
Ein paar Vorgaben für die Aufgaben:
scikit-learn oder numpy.# do not edit this cell
import warnings
warnings.filterwarnings("ignore")
Laden und visualisieren Sie den Datensatz /data/house_data.csv.
#imports
import numpy as np
import pandas as pd
import math
import matplotlib
from matplotlib import pyplot as plt
from matplotlib import cm
%matplotlib inline
import seaborn as sns
from sklearn.linear_model import LinearRegression as skLinearRegression
datapath = '/data/house_data.csv'
df = pd.read_csv(datapath)
sns.set_style('whitegrid')
_ = sns.pairplot(df)
Man sieht auf dem Pairplot auf der X sowie auf der Y-Achse (nachfolgend "grosse X & Y-Achsen genannt) alle features. Es wird als Matrix dargestellt. In den jeweiligen Kreuzungen sieht man wiederum ein Diagramm, welches auf der X-Achse den grossen X-Achsenwert und auf der Y-Achse den grossen Y-Achsenwert hat.
Schätzen Sie price mit einer linearen Regression aus dem Attribut sqft_living mit dem im Unterricht besprochenen Gradient-Descent-Verfahren. Implementieren Sie dazu eine Estimator-Klasse, analog zu scikit-learn. Visualisieren Sie die Kosten per Iteration um sicherzustellen, dass das Verfahren korrekt funktioniert.
Visualisieren Sie die Regressionsgerade zusammen mit allen Datenpunkten.
Erstellen Sie weiter einen Tukey-Ascombe-Plot und interpretieren Sie diesen.
(Detaillierte Anleitung dazu, falls nötig, finden sie hier.)
class LinearRegressionGradientDescent(object):
'''
Implements a Gradient Descent Optimization for a Simple
Linear Regression Problem.
'''
def __init__(self, init_coef=(0., 0.), epsilon=0.00001, maxsteps=1000,
linesearch='fixed', stepsize=0.001,
alpha = 0.1, beta = 0.7):
'''
Gradient Descent Optimizer for Linear Regression.
'''
self.init_coef = init_coef
self.coef_ = np.array(self.init_coef)
self._nsteps = 0
self.stepsize = stepsize
self.epsilon = epsilon
self.maxsteps = maxsteps
self.linesearch = linesearch
self.normgrad_ = []
self.cost_function_values = []
self.lsgCalculations = []
def fit(self, X, y, verbose=False, btAlpha = 0.01, btBeta = 0.38, collect_cost_function_values = True):
'''
Fits the coefficients beta of a linear regression problem
to the dataset (X, y).
'''
self.normgrad_ = []
self.cost_function_values = [];
self._nsteps = 0
while not self._do_stop(X, y):
neg_grad = -self.least_squares_gradient(self.coef_, X, y)
if self.linesearch == 'fixed':
self.coef_ += self.stepsize*neg_grad
elif self.linesearch == 'backtracking':
stepsize = self.backtracking_linesearch(X, y, alpha = btAlpha, beta = btBeta, gradFx=-neg_grad)
self.coef_ += stepsize*neg_grad
if (collect_cost_function_values):
self.cost_function_values.append(self.cost_function(X,y,self.coef_))
if verbose:
print('step {s} : neg grad = {ng}, normgrad: {nog}, coef_ = {c}, stepsize = {ss}'.format(
s=self._nsteps, ng=neg_grad, c=self.coef_,
nog=np.linalg.norm(neg_grad),
ss=self.stepsize if self.linesearch == 'fixed' else stepsize ))
self._nsteps += 1
print('Done fitting!')
print('n steps : ', self._nsteps)
print('norm gradient : ', self.normgrad_[-1])
return self
def fit_batching(self, X, y, btAlpha = 0.01, btBeta = 0.9, batchSize=10, numberOfBatches=10, printDiagramms=False, verbose=False):
'''
fittet zuerst batchsize, dann 2*batchsize, dann 3*batchsize bis numberOfBatches erreicht
und ruft dann fit(X,y,..) auf um so den gesamten Datensatz zu fitten
'''
if (verbose):
print('starting fit_batching')
totalCount = y.shape[0]
allData = X
allOutput = y
dataCount = 0
if (totalCount < numberOfBatches * batchSize):
numberOfBatches = (int)(totalCount / batchSize)
if (verbose):
print('number of batches changed! new numberOfbatches: ', numberOfBatches)
for i in range(numberOfBatches):
indezes = np.random.choice(allData.shape[0], batchSize, replace=False)
newSampleData = allData[indezes]
newSampleOutput = allOutput[indezes]
allData = np.delete(allData, indezes, 0)
allOutput = np.delete(allOutput, indezes, 0)
if (i == 0):
sampleData = newSampleData
sampleOutput = newSampleOutput
else:
sampleData = np.vstack((sampleData, newSampleData))
sampleOutput = np.append(sampleOutput, newSampleOutput)
self._nsteps = 0
while not self._do_stop(sampleData, sampleOutput):
neg_grad = -self.least_squares_gradient(self.coef_, sampleData, sampleOutput)
stepsize = self.backtracking_linesearch(sampleData, sampleOutput, alpha = btAlpha, beta = btBeta, gradFx=-neg_grad)
self.coef_ += stepsize*neg_grad
self._nsteps += 1
if (verbose):
print('done fitting ',sampleOutput.shape[0], ' data, steps: ',self._nsteps,', coefs: ', self.coef_)
print('norm gradient : ', self.normgrad_[-1])
if (printDiagramms):
b, a = self.coef_
fig, ax = plt.subplots(figsize=(5,5))
_ = ax.plot(sampleData[:,1], sampleOutput, 'bo', label='Datenpunkte')
_ = ax.plot(sampleData[:, 1], a*sampleData[:, 1]+b, '-r', label='Regressionsgerade')
self._nsteps = 0
self.fit(X, y, verbose=verbose, btAlpha=btAlpha, btBeta=btBeta)
def fit_minibatch(self, X, y, verbose=False, batchsize = 10):
'''
minibatch fitting, welche mit generationen (Epochen) arbeitet. maxstepsize wird als maxGenerations verstanden.
'''
if (verbose):
print('starting minibatch')
numberOfData = y.shape[0]
oldGradAvg = None
currGradAvg = None
generation = 0
grads = []
old_coefs = self.init_coef
self.cost_function_values.append(self.cost_function(X,y,self.coef_))
lastcf = self.cost_function_values[0]
while (generation == 0) or ((abs(self.cost_function_values[-1] - lastcf) > self.epsilon)
and (generation < self.maxsteps)):
#dont misfit
if (self.cost_function_values[-1] > lastcf):
if (verbose):
print('found a misfit')
self.coef_ = old_coefs
self.cost_function_values.append(lastcf) # for graph
old_coefs = self.coef_
permutatedIndexes = np.random.permutation(numberOfData)
if (verbose):
print('permutated Indexes: ', permutatedIndexesc)
lastcf =self.cost_function_values[-1]
for i in range(math.ceil(numberOfData / batchsize)):
batchIndex = permutatedIndexes[i*batchsize: (i+1)*batchsize]
batchX = X[batchIndex]
batchY = y[batchIndex]
neg_grad = -self.least_squares_gradient(self.coef_, batchX, batchY)
n = self.backtracking_linesearch(batchX, batchY, alpha = 0.1, beta = 0.9, gradFx=-neg_grad)
if (verbose):
print('grad: {grad} n: {n} coef: {coef}'.format(grad=neg_grad, n=n, coef=self.coef_))
self.coef_ += n * neg_grad
generation += 1
if (verbose):
print('done fitting generation',generation)
print('coefs', self.coef_)
self.cost_function_values.append(self.cost_function(X,y,self.coef_))
if (verbose):
print('last cost function value was: ', lastcf, 'curr cost function value is: ', self.cost_function_values[-1])
print('done fitting after {gen} generations'.format(gen=generation))
def _do_stop(self, X, y):
neg_grad = -self.least_squares_gradient(self.coef_, X, y)
norm_neg_grad = np.linalg.norm(neg_grad)
self.normgrad_.append(norm_neg_grad)
tiny = norm_neg_grad < self.epsilon
manysteps = self._nsteps >= self.maxsteps
return (tiny or manysteps)
def predict(self, X):
return X.dot(self.coef_)
@staticmethod
def least_squares_gradient(coef, X, y):
'''
Calculates the least squares gradient at position coef
for the dataset (X, y).
'''
grad = []
for idx in range(X.shape[1]):
grad.append(((y - coef.dot(X.T))*X[:, idx]).sum())
return -2*np.array(grad)
def score(self, X, y):
predicted = self.predict(X)
residuals = y - predicted
return 1. - ((residuals - residuals.mean())**2).sum()/((y-y.mean())**2).sum()
@staticmethod
def cost_function(X, y, beta):
return ((y - X.dot(beta))**2).sum()
def backtracking_linesearch(self, X, y, alpha = 0.01, beta = 0.38, gradFx = None):
#n = 1
#solange f(x-n*grad f(x)) > f(x) - alpha * n * grad fx)^T * Grad f(x)
# n = beta * n
coef = self.coef_
n = 1
fx = self.cost_function(X,y,coef)
if (gradFx is None):
gradFx = self.least_squares_gradient(coef,X,y);
rightPart1 = alpha * np.dot(gradFx.transpose(), gradFx)
while (self.cost_function(X,y,coef - n * gradFx) > (fx - n* rightPart1)):
n = beta * n;
return n;
def print_learningCurve(self):
fig, ax = plt.subplots(figsize=(10,6))
_ = ax.plot(self.cost_function_values, label='cost function Werte');
_ = ax.set_ylabel('cost function Werte')
_ = ax.set_xlabel('Anzahl Schritte')
_ = ax.set_title('Lern-Kurve')
_ = ax.legend()
return fig, ax
df.columns
## DATA ##
trainData=df.sample(frac=0.8, random_state=200)
testData=df.drop(trainData.index)
data_X = np.array(trainData['sqft_living'])
X_woI = np.expand_dims(data_X, axis=1) #X without Intercept
X = data_X
X = np.array([np.ones(X.shape), X]).T
y = np.array(trainData['price'])
y_log = np.log(y)
X_LAT_LONG = np.array(trainData[['lat','long']])
data_zip = np.array(trainData['zipcode'])
data_yr_built = np.array(trainData['yr_built'])
data_bedrooms = np.array(trainData['bedrooms'])
data_a8 = np.array(trainData[['bathrooms', 'grade']])
global_maxStepSize = 10**4
## TEST DATA ##
data_X_test = np.array(testData['sqft_living'])
X_woI_test = np.expand_dims(data_X_test, axis=1) #X without Intercept
X_test = data_X_test
X_test = np.array([np.ones(X_test.shape), X_test]).T
y_test = np.array(testData['price'])
y_log_test = np.log(y_test)
X_LAT_LONG_test = np.array(testData[['lat','long']])
data_zip_test = np.array(testData['zipcode'])
data_yr_built_test = np.array(testData['yr_built'])
data_bedrooms_test = np.array(testData['bedrooms'])
data_a8_test = np.array(testData[['bathrooms', 'grade']])
## STD ##
def standardize(npArr, mean = None, std = None):
if (mean is None):
mean = npArr.mean(axis=0)
if (std is None):
std = npArr.std(axis=0)
npArr_std = (npArr - mean)/std
return npArr_std, mean, std
## DATA ##
Xstd, Xmean, X_std = standardize(data_X)
Xstd = np.array([np.ones(Xstd.shape), Xstd]).T
data_a8std, data_a8mean, data_a8_std = standardize(data_a8)
## TEST DATA ##
Xstd_test, _, _ = standardize(data_X_test, Xmean, X_std)
Xstd_test = np.array([np.ones(Xstd_test.shape), Xstd_test]).T
data_a8std_test, _, _ = standardize(data_a8_test, data_a8mean, data_a8_std)
lr = LinearRegressionGradientDescent(init_coef=(0., 0.), linesearch='backtracking', epsilon=0.00001, maxsteps=global_maxStepSize)
_ = lr.fit(X, y, btAlpha=0.1, btBeta=0.7)
print('normal coefs', lr.coef_)
print('normal score', lr.score(X,y))
lr_std = LinearRegressionGradientDescent(init_coef=[0.,0.],linesearch='backtracking',epsilon=0.00001,maxsteps=global_maxStepSize)
_ = lr_std.fit(Xstd, y)
print('std coef', lr_std.coef_)
print('std score', lr_std.score(Xstd, y))
Da die Cost-function viel steiler ausfällt, wenn man standardisiert, entscheide ich mich dafür dies zu tun. Ansonsten ist die "Schale" extrem flach und selbst mit Backtracking kommt man nur sehr langsam an das Minimum.
sklr = skLinearRegression(fit_intercept=False);
sklr.fit(X,y);
print('sklr normal coefs', sklr.coef_)
print('sklr normal score', sklr.score(X, y))
sklr_std = skLinearRegression(fit_intercept=False);
sklr_std.fit(Xstd,y);
print('sklr std coefs', sklr_std.coef_)
print('sklr std score', sklr_std.score(Xstd, y))
def plot_regressionsgerade(a, b, X, y, xLabel='', yLabel='',title='Regressionsgerade',plotSKLR=False, ska = 0, skb = 0):
#a = lr_std.coef_
fig, ax = plt.subplots(figsize=(10,6))
_ = ax.plot(X[:,1], y, 'bo', label='Datenpunkte')
if (plotSKLR):
_ = ax.plot(X[:,1], ska*X[:,1]+skb, '--g', label='scikit-learn Regressionsgerade', alpha=0.5)
_ = ax.plot(X[:,1], a*X[:,1]+b, '-r', label='Regressionsgerade', alpha=0.5)
_ = ax.set_xlabel(xLabel)
_ = ax.set_ylabel(yLabel)
_ = ax.set_title(title)
_ = ax.legend()
return fig, ax
stdb, stda = lr_std.coef_
skb, ska = sklr_std.coef_
_ = plot_regressionsgerade(a= stda, b=stdb, X=Xstd, y=y, xLabel='std(sqft_living)', yLabel='price', plotSKLR=True, ska=ska,skb=skb)
Auf der X-Achse sieht man die standardisierten sqft_living, auf der Y-Achse den preis.
Meine Gerade liegt mit der von sci-kit learn übereinander. Das deutet darauf hin, dass mein Estimator für diesen Fall korrekt funktioniert.
_ = lr_std.print_learningCurve()
Auf der X-Achse sieht man die anzahl Schritte, auf der Y-Achse die Kostenfunktionswerte.
Man sieht hier, dass die Werte zuerst sehr stark abnehmen und dann nur noch langsam.
# test score with test-data
print('testdata score', lr_std.score(Xstd_test, y_test))
print('traindata score', lr_std.score(Xstd, y))
Der Score für die Test sowie für die Trainings-Daten liegen nahe beieinander, daraus lässt sich schliessen, dass der fit gut ist und kein overfitting betrieben wurde.
def print_tukey(predicted, residuals, ylim = None):
fig, ax = plt.subplots(figsize=(6,6))
_ = ax.plot(predicted, residuals, 'bo', alpha=0.2, label='residuals')
_ = ax.grid(True)
_ = ax.set_title('Tukey-Ascombe Plot')
_ = ax.set_ylabel('r')
_ = ax.set_xlabel('$\hat y$')
if (ylim is not None):
_ = ax.set_ylim(-ylim,ylim)
_ = ax.grid(b=True, which='major', color='black', linestyle='-')
_ = ax.grid(b=True, which='minor', color='r', linestyle='--')
_ = ax.legend()
return fig, ax
predicted = lr.predict(X)
residuals = y - predicted
ylim = 2000000
_ = print_tukey(predicted, residuals, ylim)
Bei dem Tukey-Ascome Plot sieht man auf der X Achse die vorhergesagten Werte des Modells. Auf der Y Achse sieht man die Abweichung vom Tatsächlichen Wert.
Man sieht, dass sich ein Trichter bildet, das deutet darauf hin, dass man logarithmisch skalieren soll.
Wiederholen Sie alle Schritte aus Teilaufgabe 2 mit der transformierten Output-Variablen log(price).
Erstellen Sie zusätzlich ein Histogramm der Residuen und interpretieren Sie dieses bezüglich Verteilung und Standardabweichung.
lr_log_std = LinearRegressionGradientDescent(init_coef=[0.,0.], linesearch='backtracking', epsilon=0.00001, maxsteps=global_maxStepSize)
_ = lr_log_std.fit(Xstd, y_log, btAlpha = 0.1, btBeta = 0.8)
stdb, stda = lr_log_std.coef_
sklr_log = skLinearRegression(fit_intercept=False)
sklr_log.fit(Xstd,y_log)
skb, ska = sklr_log.coef_
_ = plot_regressionsgerade(stda, stdb, Xstd, y_log, xLabel='std(sqft_living)', yLabel='log(price)', plotSKLR=True, ska=ska,skb=skb)
print('my coefs',lr_log_std.coef_)
print('my score training data',lr_log_std.score(Xstd, y_log))
print('my score test data',lr_log_std.score(Xstd_test, y_log_test))
Auf der X-Achse sieht man wieder den standardisierten Wert von sqft_living, auf der Y-Achse den logarithmus vom Preis.
Mein Estimator hat wieder sehr ähnlich wie sci-kit learn gefittet, was darauf schliessen lässt, dass mein Estimator korrekt funktioniert für diese Daten.
_ = lr_log_std.print_learningCurve()
Wie bereits vorher, sieht man auf der X-Achse die anzahl Schritte und auf der Y-Achse die Kostenfunktionswerte. Man sieht auch hier, dass die Werte schnell sinken und die Funktion streng monoton fallend ist.
#tukey ascombe plot
predicted = lr_log_std.predict(Xstd)
residuals = y_log- predicted
_ = print_tukey(predicted, residuals)
Auf diesem Tukey Plot sieht man, dass die Residuen viel besser verteilt sind als beim vorherigen, es ist kein einduetiger Trend mehr ausmachbar. Da es immer noch viel Abweichung hat, muss man eventuell weitere Variabeln hinzuziehen um das Resultat genauer zu bestimmen.
def draw_histogram(data, title, dataLabel='Residuen', printMeanLine=True, printStdLines=True, printMedian=True, bins = 40):
fig, ax = plt.subplots(figsize=(8, 5))
_ = ax.hist(data, bins=bins, label=dataLabel)
mean = data.mean()
if (printMeanLine):
_ = ax.axvline(mean, color='y', linewidth=2, label='Durchscnitt {mean}'.format(mean=mean))
if (printMedian):
median = np.median(data)
_ = ax.axvline(median, color='c',linewidth=2,label='Median: {median}'.format(median=median))
if (printStdLines):
std = data.std()
_ = ax.axvline(mean-std, color='r', linewidth=2, label='Durchschnitt-std')
_ = ax.axvline(mean+std, color='r', linewidth=2, label='Durchschnitt+std')
_ = ax.set_xlabel(dataLabel)
_ = ax.set_ylabel('Anzahl')
_ = ax.set_title(title)
_ = ax.legend()
return fig, ax;
residualsTrainingA3 = y_log - lr_log_std.predict(Xstd)
residualsTestA3 = y_log_test - lr_log_std.predict(Xstd_test)
_ = draw_histogram(residualsTrainingA3, 'Histogramm der Residuen Trainings Daten')
_ = draw_histogram(residualsTestA3, 'Histogramm der Residuen Test Daten')
print('Standardabweichung Residuen Trainings Daten A3', residualsTrainingA3.std())
print('Standardabweichung Residuen Test Daten A3', residualsTestA3.std())
Bei dem Histogramm sieht man die Abweichung auf der X-Achse und auf der Y-Achse, wie viele Datensätze in dem Eimer sind.
Das Diagramm ist relativ symetrisch um 0 herum. Das äussert sich auch durch einen sehr kleinen Durchschnitt (fast 0) und einen kleinen Median, der nahe beim durchschnitt liegt.
Die Standardabweichung zeigt in welchem "Rahmen" man sich bewegt, die meisten Daten liegen in dem "Schlauch" zwischen Durchschnitt + std und Durchschnitt - std
Berechnen sie den mean absolute percentage error
\begin{equation} \text{mean-ape} = \frac{1}{N}\sum_{i=1}^{N}\frac{\left| y^{(i)} - l(\beta, x^{(i)}) \right|}{\left| y^{(i)} \right|} \cdot 100 \end{equation}
für die Modelle aus Aufgabe 2 und 3.
Visualisieren sie die Verteilung des percentage error's der beiden Modelle mittels Histogrammen und diskutieren Sie diese.
def ape(y,predicted):
return (abs(y-predicted) / abs(y)) * 100
def mean_ape(y,predicted):
apes = []
for i in range(y.shape[0]):
temp = ape(y[i],predicted[i])
apes.append(temp)
return np.average(apes), apes
residuals_a2 = y-lr_std.predict(Xstd)
residuals_a2_t = y_test-lr_std.predict(Xstd_test)
residuals_a3 = y_log-lr_log_std.predict(Xstd)
residuals_a3_t = y_log_test-lr_log_std.predict(Xstd_test)
print('Residuen Durchschnitt A2',(residuals_a2).mean())
print('Residuen Durchschnitt A3',(residuals_a3).mean())
a2_MAPE, a2_apes = mean_ape(y, lr_std.predict(Xstd))
a3_MAPE, a3_apes = mean_ape(y_log, lr_log_std.predict(Xstd))
print('mape ohne log (Aufgabe 2)', a2_MAPE, '%')
print('mape log (Aufgabe 3)', a3_MAPE, '%')
a2_MAPE_t, a2_apes_t = mean_ape(y_test, lr_std.predict(Xstd_test))
a3_MAPE_t, a3_apes_t = mean_ape(y_log_test, lr_log_std.predict(Xstd_test))
print('mape ohne log (Aufgabe 2) mit Testdaten', a2_MAPE_t, '%')
print('mape log (Aufgabe 3) mit Testdaten', a3_MAPE_t, '%')
_ = draw_histogram(residuals_a2, 'Histogramm der Residuen (A2)', bins=40)
# _ = draw_histogram(residuals_a2_t, 'Histogramm der Residuen (A2) mit Testdaten')
_ = draw_histogram(residuals_a3, 'Histogramm der Residuen (A3)', bins=40)
# _ = draw_histogram(residuals_a3_t, 'Histogramm der Residuen (A3) mit Testdaten')
_ = draw_histogram(np.array(a2_apes),
'apes Aufgabe 2',
'apes',
printStdLines=False,
printMeanLine=True,
printMedian=True,
bins=100
)
# _ = draw_histogram(np.array(a2_apes_t),
# 'apes Aufgabe 2 mit Testdaten',
# 'apes',
# printStdLines=False,
# printMeanLine=True,
# printMedian=True,
# bins=100
# )
_ = draw_histogram(np.array(a3_apes),
'apes Aufgabe 3',
'apes',
printStdLines=False,
printMeanLine=True,
printMedian=True,
bins=100)
# _ = draw_histogram(np.array(a3_apes_t),
# 'apes Aufgabe 3 mit Testdaten',
# 'apes',
# printStdLines=False,
# printMeanLine=True,
# printMedian=True,
# bins=100)
Auf den ersten beiden Histogrammen sieht man auf der X-Achse die Residuen und auf der Y-Achse die Anzahl. Es werden jeweils 40 Eimer erstellt, in denen die Daten verräumt werden.
Man sieht, dass die Histogramme symetrisch verteilt sind, was darauf hindeutet, dass das Modell gut gefittet wurde. Weiter sieht man, dass der Meidan und der Durchscnitt nahe beieinander liegen, das deutet weiter darauf hin, dass die Verteilung symetrisch ist.
Bei Aufgabe 3 sieht man, dass der Durchschnitt näher bei 0 liegt als bei Aufgabe 2.
Bei den zweiten beiden Histogrammen sieht man auf der X-Achse die Percentage-Errors und auf der Y-Achse die Anzahl. Es werden wieder jeweils 100 Buckets erstellt, in denen die Daten verräumt werden.
Bei Aufgabe 2 sieht man, dass die Daten langsam abnehmen, der Durschnitt jedoch bei 36% liegt. der Median ist kleiner als der Durchscnitt, was auf ein paar Ausreisser hindeutet, welche den Durchschnitt erhöhen.
Bei Aufgabe 3 sieht man einen steileren Abfall, Durchschnitt und Median liegen hier näher beieinander. Dies deutet darauf hin, dass die Ausreisser von Aufgabe 2 durch das logarithmische Transformieren ebenfalls gut gefittet werden konnten.
Visualisieren Sie mit einer geeigneten Farbskalierung die Output-Variablen log(‘price‘) als Punkte in den geografischen Koordinaten 'long'/'lat' (geografische Länge/Breite) und diskutieren Sie das Ergebnis.
fig, ax = plt.subplots(figsize=(10, 8))
sc = ax.scatter(X_LAT_LONG[:,1],X_LAT_LONG[:,0], c=y_log, cmap=cm.Reds)
ax.set_xlabel('lat')
ax.set_ylabel('long')
ax.set_title('Preis auf der Karte')
_ = fig.colorbar(sc)
Wiederholen Sie Aufgabe 5 mit dem Attribut 'zipcode' anstelle von log('price') und interpretieren Sie den Einfluss von 'zipcode' auf den Immobilienpreis.
fig, ax = plt.subplots(figsize=(8, 8))
sc = ax.scatter(X_LAT_LONG[:,1],X_LAT_LONG[:,0], c=data_zip, cmap=cm.Reds)
ax.set_xlabel('lat')
ax.set_ylabel('long')
ax.set_title('Zipcode auf der Karte 1')
_ = fig.colorbar(sc)
fig, ax = plt.subplots(figsize=(8, 8))
sc = ax.scatter(X_LAT_LONG[:,1],X_LAT_LONG[:,0], c=data_zip, cmap=cm.tab20c)
ax.set_xlabel('lat')
ax.set_ylabel('long')
ax.set_title('Zipcode auf der Karte 2')
_ = fig.colorbar(sc)
Auf dem ersten Plot sieht man die Postleizahlen auf der Karte. Höhere Postleizahlen haben eine rötere Farbe als tiefere. Wenn man diesen Plot mit dem Preis-Plot vergleicht, sieht man keine zusammenhänge. Also eine höhere Postleizahl bedeutet nicht einen höheren Preis. (Dies würde man auch auf dem ersten (seaborn pairplot) plot erkennen und ist hier nur vollständigkeitshalber aufgeführt)
Auf dem zweiten Plot sieht man die Gebiete besser. Wenn man den Plot aus Aufgabe 5 mit diesem vergleicht, sieht man, dass gewisse Postleizahlen eher teuerer sind als andere. Daher empfielt es sich, den Zip-code ebenfalls als Feature hinzuzuziehen um genauere Schätzungen zu machen.
Feature Engineering: Codieren Sie das kategorische Attribut 'zipcode' mittels One-Hot Encoding und schätzen Sie log('price') mit einer neuen linearen Regression aus 'sqft_living' und den neu kreierten Attributen.
Erstellen Sie ein Histogramm der Residuen und diskutieren Sie dieses.
Wie verhält sich der Mittelwert, bzw. die Standardabweichung im Vergleich zu Aufgabe 3?
X_ohe = Xstd;
coef_init = lr_log_std.coef_ # reuse old coef, faster fitting :)
for zip in np.unique(data_zip):
encoded = data_zip == zip
encoded = np.array(encoded).astype('int')
encoded = np.expand_dims(encoded, axis=1)
X_ohe = np.hstack((X_ohe, encoded))
coef_init = np.hstack((coef_init, 0))
## TEST DATA ##
X_ohe_test = Xstd_test
for zip in np.unique(data_zip_test):
encoded = data_zip_test == zip
encoded = np.array(encoded).astype('int')
encoded = np.expand_dims(encoded, axis=1)
X_ohe_test = np.hstack((X_ohe_test, encoded))
lr_minibatch_a7 = LinearRegressionGradientDescent(init_coef=coef_init,linesearch='backtracking',epsilon=10**-1,maxsteps=100)
_ = lr_minibatch_a7.fit_minibatch(X_ohe, y_log, batchsize=2000, verbose=False)
print('minibatch coef', lr_minibatch_a7.coef_)
print('minibatch score', lr_minibatch_a7.score(X_ohe, y_log))
fig, ax = lr_minibatch_a7.print_learningCurve()
ax.set_xlabel('# Generationen + mismatchkorrekturen')
ax.set_title('Lern-Kurve mini-batching A7')
lr_after_minibatch_a7 = LinearRegressionGradientDescent(init_coef=lr_minibatch_a7.coef_, linesearch='backtracking', epsilon=10**-2, maxsteps=global_maxStepSize)
_ = lr_after_minibatch_a7.fit(X_ohe, y_log)
print('after minibatch coef', lr_after_minibatch_a7.coef_)
print('after minibatch score', lr_after_minibatch_a7.score(X_ohe, y_log))
fig, ax = lr_after_minibatch_a7.print_learningCurve()
Auf der ersten Lern-Kurve sieht man das mini-batching. Im Gegensatz zu der gewöhnlichen Lern-Kurve kann diese wieder zunehmen. Allerdings verwerfe ich in der Funktion diese Zunahmen wieder und sie werden nur auf dem plot dargestellt. Es werden immer die Koeffizienten genommen, welche die minimale Kostenfunktion gegeben haben.
Da die Permutationen zufällig gewählt werden, kann es sein, dass bei Ihrer Ausführung der Graph etwas anders aussieht. Ebenfalls kann es sein, dass es erst nach 100 Generationen (max) terminiert.
first_cfv = lr_minibatch_a7.cost_function_values
second_cfv = lr_after_minibatch_a7.cost_function_values
both_cost_values = []
for cfv in first_cfv:
both_cost_values.append(cfv)
for cfv in second_cfv:
both_cost_values.append(cfv)
#print(first_cfv)
fig, ax = plt.subplots(figsize=(10,6))
_ = ax.plot(both_cost_values, label='cost function Werte');
_ = ax.set_ylabel('cost function Werte')
_ = ax.set_xlabel('Anzahl Generationen/Schritt/')
_ = ax.set_title('Lern-Kurve')
_ = ax.legend()
Dies sind die beiden Lernkurven aneinandergehängt. Man sieht hier, dass ein grosser Teil der Arbeit durch das mini-batching erledigt werden konnte und das normale fitten nur noch eine Feinpolitur war.
#old code
#lr_Zip = LinearRegressionGradientDescent(init_coef=coef_init, linesearch='backtracking', epsilon=0.00001, maxsteps=global_maxStepSize)
#_ = lr_Zip.fit_batching(X_ohe, y_logstd, batchSize=10, numberOfBatches=15, verbose=True)
print('score', lr_after_minibatch_a7.score(X_ohe, y_log))
residualsTrainingA7 = y_log - lr_after_minibatch_a7.predict(X_ohe)
residualsTestA7 = y_log_test - lr_after_minibatch_a7.predict(X_ohe_test)
_ = draw_histogram(data=residualsTrainingA3, title='Residuen A3 (Trainings-daten)')
_ = draw_histogram(data=residualsTrainingA7, title='Residuen A7 (Trainings-daten)')
_ = draw_histogram(data=residualsTestA7, title='Residuen A7 (Test-daten)')
print('Standardabweichung Residuen Trainings Daten A3', residualsTrainingA3.std())
print('Standardabweichung Residuen Test Daten A3', residualsTestA3.std())
print('Standardabweichung Residuen Trainings Daten A7', residualsTrainingA7.std())
print('Standardabweichung Residuen Test Daten A7', residualsTestA7.std())
fig, ax = print_tukey(lr_log_std.predict(Xstd), residualsTrainingA3)
_ = ax.set_title('Tukey-Ascombe Plot A3')
fig, ax = print_tukey(lr_after_minibatch_a7.predict(X_ohe), residualsTrainingA7)
_ = ax.set_title('Tukey-Ascombe Plot A7')
Auf dem Histogramm sieht man in der X-Achse die Residuen und auf der Y-Achse die Anzahl der Einträge im Eimer. Man sieht hier, dass der Durchschnitt wieder bei praktisch 0 liegt und der Median sehr nahe bei 0. Die Standardabweichung konnte im Gegensatz zu Aufgabe 3 fast halbiert werden, das Modell streut nun nicht mehr so stark. Das lässt sich auch durch einen Tukey-Ascombe Plot visualisieren, hier ist der Bereich, in dem sich die Residuen bewegen von ugnefähr +-1 auf ugnefähr +-0.5 gesunken
Der Durchschnitt liegt dieses mal im positiven bereich, vorher war dieser negativ. Das Bedeutet, dass das Modell die Preise im Durchschnitt zu hoch schätzt, und vorher zu tief.
## checking
sklr = skLinearRegression(fit_intercept=False);
sklr.fit(X_ohe,y_log);
print('score sklr - my score: ', sklr.score(X_ohe,y_log) - lr_after_minibatch_a7.score(X_ohe, y_log))
print('sklr.coef_ - my.coef_:')
print(sklr.coef_- lr_after_minibatch_a7.coef_)
Wiederholen Sie die Aufgabe 7 mit den vier zusätzlichen Attributen ['bedrooms', 'bathrooms', 'grade', 'yr_built'] und diskutieren Sie den Unterschied zu Aufgabe 7.
#entscheidungshilfe
a8_evaluation_df = df[['bedrooms', 'bathrooms', 'grade','yr_built','price']]
sns.set_style('whitegrid')
_ = sns.pairplot(a8_evaluation_df)
Aus dem Pairplot ist ersichtlich, dass bathrooms und grade einigermassen linear abhängig sind. Bei yr_built ist keine lineare abhängigkeit zu sehen, ich werde dieses feature encoden. bei bedrooms scheint es so zu sein, als ob es eine "super" anzahl gäbe, bei der der Preis hoch ist, dieser danach jedoch wieder abfällt (evtl Hotels?). um dies zu berücksichtigen werde ich auch dieses feature encoden.
X_ohe2 = X_ohe;
coef_init = lr_after_minibatch_a7.coef_ #again, reuse coefs for faster fitting
# encode yr_built with ohe
for yr_built in np.unique(data_yr_built):
encoded = data_yr_built == yr_built
encoded = np.array(encoded).astype('int')
encoded = np.expand_dims(encoded, axis=1)
X_ohe2 = np.hstack((X_ohe2, encoded))
coef_init = np.hstack((coef_init, 0))
#encode bedrooms with ohe
for bedrooms in np.unique(data_bedrooms):
encoded = data_bedrooms == bedrooms
encoded = np.array(encoded).astype('int')
encoded = np.expand_dims(encoded, axis=1)
X_ohe2 = np.hstack((X_ohe2, encoded))
coef_init = np.hstack((coef_init, 0))
X_ohe2 = np.hstack((X_ohe2,data_a8std))
coef_init = np.hstack((coef_init, np.zeros(data_a8std.shape[1])))
#test-data
## TEST DATA ##
X_ohe2_test = X_ohe_test
# encode yr_built with ohe
for yr_built in np.unique(data_yr_built_test):
encoded = data_yr_built_test == yr_built
encoded = np.array(encoded).astype('int')
encoded = np.expand_dims(encoded, axis=1)
X_ohe2_test = np.hstack((X_ohe2_test, encoded))
#encode bedrooms with ohe
for bedrooms in np.unique(data_bedrooms_test):
encoded = data_bedrooms_test == bedrooms
encoded = np.array(encoded).astype('int')
encoded = np.expand_dims(encoded, axis=1)
X_ohe2_test = np.hstack((X_ohe2_test, encoded))
X_ohe2_test = np.hstack((X_ohe2_test,data_a8std_test))
lr_minibatch_a8 = LinearRegressionGradientDescent(init_coef=coef_init,linesearch='backtracking',epsilon=10**-1,maxsteps=100)
_ = lr_minibatch_a8.fit_minibatch(X_ohe2, y_log, batchsize=2000, verbose=False)
lr_after_minibatch_a8 = LinearRegressionGradientDescent(init_coef=lr_minibatch_a8.coef_, linesearch='backtracking', epsilon=3, maxsteps=global_maxStepSize)
_ = lr_after_minibatch_a8.fit(X_ohe2, y_log)
# print('minibatch coef', lr_minibatch_a8.coef_)
print('minibatch score', lr_minibatch_a8.score(X_ohe2, y_log))
fig, ax = lr_minibatch_a8.print_learningCurve()
ax.set_xlabel('# Generationen (+mitsmatch korrekturen)')
ax.set_title('Lern-Kurve mini-batching A8')
# print('after minibatch coef', lr_after_minibatch_a8.coef_)
print('after minibatch score', lr_after_minibatch_a8.score(X_ohe2, y_log))
print('after minibatch score Testdaten', lr_after_minibatch_a8.score(X_ohe2_test, y_log_test))
fig, ax = lr_after_minibatch_a8.print_learningCurve()
first_cfv = lr_minibatch_a8.cost_function_values
second_cfv = lr_after_minibatch_a8.cost_function_values
both_cost_values = []
for cfv in first_cfv:
both_cost_values.append(cfv)
for cfv in second_cfv:
both_cost_values.append(cfv)
#print(first_cfv)
fig, ax = plt.subplots(figsize=(10,6))
_ = ax.plot(both_cost_values, label='cost function Werte');
_ = ax.set_ylabel('cost function Werte')
_ = ax.set_xlabel('Anzahl Generationen/Schritte')
_ = ax.set_title('Lern-Kurve')
_ = ax.legend()
print(lr_after_minibatch_a8.score(X_ohe2, y_log))
residualsTrainingA8 = y_log - lr_after_minibatch_a8.predict(X_ohe2)
residualsTestsA8 = y_log_test - lr_after_minibatch_a8.predict(X_ohe2_test)
draw_histogram(data=residualsTrainingA8, title='Residuen')
fig, ax = print_tukey(lr_after_minibatch_a8.predict(X_ohe2), residualsTrainingA8)
_ = ax.set_title('Tukey-Ascombe Plot A8')
print('Standardabweichung Residuen Trainings Daten A3', residualsTrainingA3.std())
print('Standardabweichung Residuen Trainings Daten A7', residualsTrainingA7.std())
print('Standardabweichung Residuen Trainings Daten A8', residualsTrainingA8.std())
print('Standardabweichung Residuen Test Daten A8', residualsTestsA8.std())
Der Score wurde nochmals etwas besser, die Standardabweichung ebenfalls. Es wurde noch kein overfitting betrieben, die Test-daten scores & standardabweichungen stimmen immernoch mit den Trainingsdaten überein.
Im Unterschied zu A7 werden die Preise im Durchschnitt wieder niedriger geschätzt als tatsächlich.
Was müssten Sie anpassen wenn Sie den MAPE als Kostenfunktion für Gradient-Descent verwenden möchten? Was könnte passieren, wenn Sie den RSS als Kostenfunktion verwenden, jedoch am Ende das Modell mittels MAPE evaluieren?
Wenn man MAPE als Kostenfunktion verwenden würde, dann könnte man das Problem nicht mehr analytisch lösen (Wegen der Betragsfunktion).
Anzupassen wären die Funktionen least_squares_gradient und cost_function. bei RSS werden weit entfernte Datenpunkte stärker bestraft (quadratischer Abstand) als bei MAPE. Dies führt dazu, dass modelle welche z.T. starke ausschläge haben mit RSS gefittet werden, jedoch immernoch einen grossen MAPE haben, da die "Ausnahmefälle" diesen